#include <stdio.h>
typedef struct
	{
		int fila;
		int columna;
		double * data;
	} str_matrix;

/* CREA UNA MATRIZ DE mxn */
str_matrix create_matrix(int m,int n){
	str_matrix mat;
	mat.fila=m;mat.columna=n;
	mat.data=(double*)malloc(m*n*sizeof(double));
	return mat;
}
/* SETEA LOS VALORES DE LA MATRIZ CON AQUELLOS PASADOS EN EL ARRAY DE doubles */
void set_values(str_matrix * mat,double * datos){
	int cont;
	int cant;
	double * p;
	p=mat->data;
	cont=0;
	cant=mat->columna*mat->fila;
	while(cont <cant){
		*p=*datos;
		p++;
		datos++;
		cont++;
	}
}
/* SETEA TODOS LOS VALORES DE LA MATRIZ A x */
void set_to_x(str_matrix * mat,double x){
	double * puntero=mat->data;
	int aux=mat->fila*mat->columna;
	
	while (aux>0){
		*puntero=x;
		aux--;
		puntero++;
	}
}
/* SETEA TODOS LOS VALORES DE LA MATRIZ A CERO */
str_matrix zeros(int m,int n){
	str_matrix mat=create_matrix(m,n);
	set_to_x(&mat,0.0);
	return mat;
}
/* MUESTRA EN CONSOLA LA MATRIZ */
void mostrar_matrix(str_matrix * mat){
	double * puntero;
	int n;
	int m;
	puntero=mat->data;
	
	m=mat->fila;
	n=mat->columna;
	
	while (m>0){
		while(n>0){
			printf("%5f\t",*puntero);
			puntero++;
			n--;
		}
		n=mat->columna;
		m--;
		printf("\n");
	}
	
}
/* MUESTRA EN CONSOLA LA MATRIZ */
void mostrar_matrix_2(str_matrix * mat){
	double * puntero;
	int n;
	int m;
	puntero=mat->data;	
	m=mat->fila;
	n=mat->columna;
	
	while (m>0){
		while(n>0){
			printf("%5f\t",*puntero);
			puntero++;
			n--;
		}
		n=mat->columna;
		m--;
	}
	
}
/* LIBERA LA MEMORIA UTILIZADA POR EL ARRAY DINAMICO DE doubles QUE TIENE LA STRUCT str_matrix */
void destroy_matrix(str_matrix * mat){
	free(mat->data);
}
/* SUMA MATRICIAL A=A+B */
void suma_2(str_matrix * m1,str_matrix * m2){
	double * puntero1;
	double * puntero2;
	int cant;
	
	puntero1=m1->data;
	puntero2=m2->data;
	
	if(m1->columna==m2->columna && m1->fila==m2->fila){
		cant=m1->columna*m1->fila;
		
		while (cant>0){
			*puntero1=*puntero1+*puntero2;	
			puntero2++;
			puntero1++;
			cant--;	
		}
	}
	else{
		printf("suma_2: Las matrices a sumar no son de las mismas dimensiones\n");
	}
}
/* SUMA MATRICIAL RET=A+B -- NO UTILIZAR PARA CALCULAR A=A-B PORQUE SE PRODUCE PERDIDA DE PUNTEROS Y DESPERDICIO DE MEMORIA */
str_matrix suma(str_matrix * m1,str_matrix * m2){
	double * puntero1;
	double * puntero2;
	double * puntero3;
	int cant;
	str_matrix mat;
	
	puntero1=m1->data;
	puntero2=m2->data;
	
	if(m1->columna==m2->columna && m1->fila==m2->fila){
		mat=create_matrix(m1->fila,m1->columna);
		puntero3=mat.data;
		cant=m1->columna*m1->fila;
		
		while (cant>0){
			*puntero3=*puntero2+*puntero1;	
			puntero3++;
			puntero2++;
			puntero1++;
			cant--;	
		}
	}
	else{
		printf("suma: Las matrices a sumar no son de las mismas dimensiones\n");
	}
	
	return mat;
}
/* RESTA MATRICIAL A=A-B */
str_matrix resta(str_matrix * m1,str_matrix * m2){
	double * puntero1;
	double * puntero2;
	double * puntero3;
	int cant;
	str_matrix mat;
	
	puntero1=m1->data;
	puntero2=m2->data;
	
	if(m1->columna==m2->columna && m1->fila==m2->fila){
		mat=create_matrix(m1->fila,m1->columna);
		puntero3=mat.data;
		cant=m1->columna*m1->fila;
		
		while (cant>0){
			*puntero3=*puntero1-*puntero2;	
			puntero3++;
			puntero2++;
			puntero1++;
			cant--;	
		}
	}
	else{
		printf("resta: Las matrices a restar no son de las mismas dimensiones\n");
	}
	
	return mat;
}
/* RESTA MATRICIAL RET=A-B -- NO UTILIZAR PARA CALCULAR A=A-B PORQUE SE PRODUCE PERDIDA DE PUNTEROS Y DESPERDICIO DE MEMORIA */
void resta_2(str_matrix * m1,str_matrix * m2){
	double * puntero1;
	double * puntero2;
	int cant;
	
	puntero1=m1->data;
	puntero2=m2->data;
	
	if(m1->columna==m2->columna && m1->fila==m2->fila){
		cant=m1->columna*m1->fila;
		
		while (cant>0){
			*puntero1=*puntero1-*puntero2;	
			puntero2++;
			puntero1++;
			cant--;	
		}
	}
	else{
		printf("resta_2: Las matrices a restar no son de las mismas dimensiones\n");
	}
}
/* DEVUELVE LA iesima FILA DE A COMO MATRIZ DE 1xn*/
str_matrix get_fila(str_matrix * mat, int indice){
	str_matrix matrix_fila;
	int i;
	double * p1;
	double * p2;
	
	i=0;
	matrix_fila.fila=1;
	matrix_fila.columna=mat->columna;
	
	p1=(double*)malloc(mat->columna*sizeof(double));
	matrix_fila.data=p1;
	p2=mat->data;
	
	p2=p2+ (indice)*(mat->columna);
	
	while(i<mat->columna){
		*p1=*p2;
		p1++;
		p2++;
		i++;
	}
	return matrix_fila;
}
/* SUSTITUYE LA iesima FILA DE A POR LA MATRIZ FILA */
void set_fila(str_matrix * mat,str_matrix * fila,int indice){
	double * pmat;
	double * pfila;
	int i;
	pmat=mat->data;
	pfila=fila->data;
	pmat=pmat+(mat->columna*indice);
	i=0;
	while(i<mat->columna){
		*pmat=*pfila;
		pmat++;
		pfila++;
		i++;
	}	
}	
/* VALOR ABSOLUTO x=|x| */
double valor_absoluto(double d){
	if (d<0){
		d=d*-1;
	}	
	return d;
}
/* AUXILIAR PARA QUITAR CEROS DE LUGARES CRITICOS QUE GENEREN DIVISION POR CERO */
void combinar_fila(str_matrix * mat,str_matrix * inv, int indice){
	int m;
	int contador;
	double min_valor;
	short flag;
	short exito;
	str_matrix m_fila;
	str_matrix m_fila_aux;
	str_matrix inv_fila;
	str_matrix inv_fila_aux;
	
	m=mat->fila;
	flag=1;
	exito=-1;
	contador=0;
	
	while(flag==1 && contador<m){
		if(contador!=indice){
			m_fila=get_fila(mat,indice);
			m_fila_aux=get_fila(mat,contador);
			inv_fila=get_fila(inv,indice);
			inv_fila_aux=get_fila(inv,contador);			
			
			suma_2(&m_fila,&m_fila_aux);
			suma_2(&inv_fila,&inv_fila_aux);
			
			min_valor=valor_absoluto(*(m_fila.data + indice));
			
			if(min_valor>0.00001){
				flag=-1;
				exito=1;
				set_fila(mat,&m_fila,indice);
				set_fila(inv,&inv_fila,indice);
			}
			destroy_matrix(&m_fila);
			destroy_matrix(&m_fila_aux);
			destroy_matrix(&inv_fila);
			destroy_matrix(&inv_fila_aux);
		}
		contador++;
	}
	if(exito==-1){
		printf("combinar_fila: No pudo eliminarse los ceros de la diagonal. Error de division por cero inminente.\n");
	}
}
/* AUXILIAR PARA QUITAR CEROS DE LUGARES CRITICOS QUE GENEREN DIVISION POR CERO */
void combinar_fila_2(str_matrix * mat,str_matrix * inv,int y,int x){
	int m;
	int contador;
	double min_valor;
	short flag;
	short exito;
	str_matrix m_fila;
	str_matrix m_fila_aux;
	str_matrix inv_fila;
	str_matrix inv_fila_aux;
	
	m=mat->fila;
	flag=1;
	exito=-1;
	contador=0;
	
	while(flag==1 && contador<m){
		if(contador!=y){
			m_fila=get_fila(mat,y);
			inv_fila=get_fila(inv,y);
			m_fila_aux=get_fila(mat,contador);
			inv_fila_aux=get_fila(inv,contador);
			suma_2(&m_fila,&m_fila_aux);
			suma_2(&inv_fila,&inv_fila_aux);
			
			min_valor=valor_absoluto(*(m_fila.data + x));
			if(min_valor>0.00001){
				flag=-1;
				exito=1;
				set_fila(mat,&m_fila,y);
				set_fila(inv,&inv_fila,y);
			}
			destroy_matrix(&m_fila);
			destroy_matrix(&m_fila_aux);
			
			destroy_matrix(&inv_fila);
			destroy_matrix(&inv_fila_aux);
		}
		contador++;
	}
	if(exito==-1){
		printf("combinar_fila_2: No pudieron eliminarse los ceros. Error de division por cero inminente.\n");
	}
}
/* PRODUCTO RET = M1*M2 -- NO UTILIZAR PARA HACER M1=M1*M2 PORQUE SE PRODUCE PERDIDA DE PUNTEROS Y DESPERDICIO DE MEMORIA*/
str_matrix prod(str_matrix * m1,str_matrix * m2){
	int m;
	int n;
	int p;
	int q;
	double * p1;
	double * p2;
	double * p3;
	int aux_f;
	int aux_c;
	int aux_c2;
	str_matrix mat;
	
	aux_f=0;
	aux_c=0;
	aux_c2=0;

	m=m1->fila;n=m1->columna;
	p=m2->fila;q=m2->columna;
	p1=m1->data;
	p2=m2->data;

	if(n==p){
		mat=zeros(m,q);
		p3=mat.data;
		
		while (aux_f<m){
			p1=m1->data;
			p1=p1+(aux_f*n);
			p2=m2->data;

			while(aux_c2<q){

				p2=p2+(aux_c2);
				
				while (aux_c<p){
				
					*p3=*p3+ ((*p1)*(*p2));
					p2=p2+q;
					p1++;
					aux_c++;
				}
				p1=m1->data;
				p1=p1+(aux_f*n);
				p2=m2->data;
				aux_c=0;
				aux_c2++;
				p3++;
			}
			aux_c2=0;
			aux_f++;
		}
	}
	else{
		printf("prod: Matrices no compatibles\n");
	}
	
	return mat;
}
/* PRODUCTO M1= M1*M2 */
void prod_2(str_matrix * m1,str_matrix * m2){
	int m;
	int n;
	int p;
	int q;
	double * p1;
	double * p2;
	double * p3;
	int aux_f;
	int aux_c;
	int aux_c2;
	str_matrix mat;
	
	aux_f=0;
	aux_c=0;
	aux_c2=0;
	
	m=m1->fila;n=m1->columna;
	p=m2->fila;q=m2->columna;
	p1=m1->data;
	p2=m2->data;
	
	if(n==p){
		mat=zeros(m,q);
		p3=mat.data;
		
		while (aux_f<m){
			p1=m1->data;
			p1=p1+(aux_f*n);
			p2=m2->data;
			
			while(aux_c2<q){
				
				p2=p2+(aux_c2);
				
				while (aux_c<p){
					
					*p3=*p3+ ((*p1)*(*p2));
					p2=p2+q;
					p1++;
					aux_c++;
				}
				p1=m1->data;
				p1=p1+(aux_f*n);
				p2=m2->data;
				aux_c=0;
				aux_c2++;
				p3++;
			}
			aux_c2=0;
			aux_f++;
		}
	}
	else{
		printf("prod_2: Matrices no compatibles\n");
	}
	destroy_matrix(m1);
	*m1=mat;
}
/* PRODUCTO M2= M1*M2 */
void prod_3(str_matrix * m1,str_matrix * m2){
	int m;
	int n;
	int p;
	int q;
	double * p1;
	double * p2;
	double * p3;
	int aux_f;
	int aux_c;
	int aux_c2;
	str_matrix mat;
	
	aux_f=0;
	aux_c=0;
	aux_c2=0;
	
	m=m1->fila;n=m1->columna;
	p=m2->fila;q=m2->columna;
	p1=m1->data;
	p2=m2->data;
	
	if(n==p){
		mat=zeros(m,q);
		p3=mat.data;
		
		while (aux_f<m){
			p1=m1->data;
			p1=p1+(aux_f*n);
			p2=m2->data;
			
			while(aux_c2<q){
				
				p2=p2+(aux_c2);
				
				while (aux_c<p){
					
					*p3=*p3+ ((*p1)*(*p2));
					p2=p2+q;
					p1++;
					aux_c++;
				}
				p1=m1->data;
				p1=p1+(aux_f*n);
				p2=m2->data;
				aux_c=0;
				aux_c2++;
				p3++;
			}
			aux_c2=0;
			aux_f++;
		}
	}
	else{
		printf("prod_3: Matrices no compatibles\n");
	}
	destroy_matrix(m2);
	*m2=mat;
}
/* PRODUCTO ESCALAR DONDE A=e*A */
void prod_escalar_2(str_matrix * mat,double num2){
	int i;
	int cant;
	double num;
	double * p;
	
	num=num2;
	cant=mat->columna*mat->fila;
	p=mat->data;
	i=0;
	while(i<cant){
		*p= (*p) * num;
		p++;
		i++;
	}
}
/* PRODUCTO ESCALAR DONDE ret=e*A */
str_matrix prod_escalar(str_matrix * mat,double num){
	str_matrix ret;
	int i;
	int cant;
	double * p1;
	double * p2;
	
	ret=create_matrix(mat->fila,mat->columna);
	cant=mat->columna*mat->fila;
	p1=mat->data;
	p2=ret.data;
	i=0;
	while(i<cant){
		*p2= (*p1) * num;
		p1++;
		p2++;
		i++;
	}
	return ret;
}
/* DIVISION ESCALAR DONDE A=A/e */
void div_escalar(str_matrix * mat,double num){
	int i;
	int cant;
	double * p;
	
	cant=mat->columna*mat->fila;
	p=mat->data;
	i=0;
	while(i<cant){
		*p= (*p) / num;
		p++;
		i++;
	}
}
/* GENERA MATRIZ IDENTIDAD DE nxn */
str_matrix eye(int n){
	str_matrix mat;
	int i;
	int j;
	double * puntero;
	
	i=0;
	j=0;
	mat=zeros(n,n);
	puntero=mat.data;
	
	while (i<n){
		while(j<n){
			if(i==j){
				*puntero=1;
			}
			puntero++;
			j++;
		}
		j=0;
		i++;
	}
	return mat;
}
/* COPIA UNA MATRIZ */
str_matrix copy_matrix(str_matrix * origen){
	int m_;
	int n_;
	int cant;
	int cont;
	str_matrix destino;
	double * p1;
	double * p2;
	
	m_=origen->fila;
	n_=origen->columna;
	cant=m_*n_;
	cont=0;
	destino=create_matrix(m_,n_);
	p1=origen->data;
	p2=destino.data;
	
	while(cont<cant){
		*p2=*p1;
		p1++;
		p2++;
		cont++;
	}
	return destino;
}
/* INVERSA DE UNA MATRIZ */
str_matrix inversa(str_matrix * original){
	int m;
	int n;
	int x;
	int y;
	int z;
	int contador;
	double num_aux;
	double zero_check;
	str_matrix m1;
	str_matrix inv;
	str_matrix aux1;
	str_matrix aux2;
	str_matrix aux3;
	str_matrix aux4;
	
	m1=copy_matrix(original);
	m=m1.fila;
	n=m1.columna;
	x=0;
	y=0;
	z=1;

	inv=eye(n);

	contador=0;
	while(contador<m){
		zero_check= *(m1.data + contador*m + contador);
		zero_check=valor_absoluto(zero_check);
		if(zero_check<0.00001){
			combinar_fila(&m1,&inv,contador);
		}
		contador++;
	}
	
	if(m==n){

		
		while(x<n){
			y=z;
			aux3=get_fila(&m1,z-1);
			aux4=get_fila(&inv,z-1);
	
			while(y<m){
				aux1=get_fila(&m1,y);
				
				zero_check=(*(aux1.data+x));
				zero_check=valor_absoluto(zero_check);
				
				if(zero_check>0.00001){
					
					aux2=get_fila(&inv,y);
					num_aux=(*(aux3.data+x))/(*(aux1.data+x));
					
					prod_escalar(&aux1,num_aux);
					prod_escalar(&aux2,num_aux);
					
					resta_2(&aux1,&aux3);
					resta_2(&aux2,&aux4);
					
					set_fila(&m1,&aux1,y);
					set_fila(&inv,&aux2,y);
					
					destroy_matrix(&aux2);
				}
				destroy_matrix(&aux1);
				y++;
			}
			x++;
			z++;
			destroy_matrix(&aux3);
			destroy_matrix(&aux4);
		}

		z=m-2,
		x=n-1;
		
		while(x>0){
			y=z;
			aux3=get_fila(&m1,z+1);
			aux4=get_fila(&inv,z+1);
			while(y>-1){

				aux1=get_fila(&m1,y);
				zero_check=(*(aux1.data+x));
				zero_check=valor_absoluto(zero_check);
				
				if(zero_check>0.00001){
					aux2=get_fila(&inv,y);
					num_aux=(*(aux3.data+x))/(*(aux1.data+x));
					
					prod_escalar(&aux1,num_aux);
					prod_escalar(&aux2,num_aux);
					
					resta_2(&aux1,&aux3);
					resta_2(&aux2,&aux4);
					
					set_fila(&m1,&aux1,y);
					set_fila(&inv,&aux2,y);
					
					destroy_matrix(&aux2);
				}
				destroy_matrix(&aux1);
				y--;
			}
			x--;
			z--;
			destroy_matrix(&aux3);
			destroy_matrix(&aux4);
		}

		x=0;
		while(x<n){
			aux1=get_fila(&m1,x);
			aux2=get_fila(&inv,x);
			div_escalar(&aux1,*(m1.data +x*n +x));
			div_escalar(&aux2,*(m1.data +x*n +x));
			set_fila(&m1,&aux1,x);
			set_fila(&inv,&aux2,x);
			
			destroy_matrix(&aux1);
			destroy_matrix(&aux2);
			
			x++;
		}
	}
	else{
		printf("inversa: La matriz no es cuadrada");
	}
	return inv;
}
/* TRANSPUESTA DE UNA MATRIZ*/
str_matrix transpuesta(str_matrix * original){
	int m;
	int n;
	int y;
	int x;

	str_matrix transp;
	double * p1;
	double * p2;
	
	m=original->fila;
	n=original->columna;

	x=0;
	y=0;
	transp=create_matrix(n,m);
	
	p1=original->data;
	
	while(y<m){
		x=0;
		while(x<n){
			p2=(transp.data + x*m + y);
			*p2=*p1;
			p1++;
			x++;
		}
		y++;
	}
	return transp;
}
/**/
